npj Systems Biology and Applications

Glycosphingolipids (GSLs) are essential components of biological membranes with important roles in cell signalling. Disrupted GSL metabolism is associated with malignancy across a range of cancers, with different GSLs implicated in distinct tumours. GSLs have potential mechanistic roles in cancer; however, their functions in Lower Grade Gliomas (LGGs) remain poorly understood. We present ensemble machine learning approaches using transcriptomic data from LGG, combined with GSL-specific metabolic simulations, to predict survival outcomes. The ensemble approach demonstrates effective risk stratification for LGG patients based on GSL gene expression. Pathway analysis of model-derived risk groups highlighted potential association of GSLs with cell motility, division and Wnt signalling in LGG pathology. Given the strong performance of machine learning approaches to predict outcomes and that GSLs are shed into the tumour microenvironment, GSL-based diagnostics and prognostics may prove clinically beneficial. A Python package enabling GSL-specific metabolic modelling and risk prediction from RNA-seq data is provided.

Anticancer drug resistance is challenging to overcome because it can arise through both intrinsic and acquired mechanisms, each driven by distinct cellular machinery. In particular, there is a sharp need for therapies that target hormone-insensitive prostate tumors due to the growing incidence of castration-resistant prostate cancer. Optimizing the pathways that regulate apoptosis in prostate cancer offers a promising strategy to induce apoptosis and inhibit tumor progression, since these mechanisms do not depend on hormonal signaling. Here, we identified strategies to enhance apoptosis in prostate cancer cells. We used several computational tools (including sensitivity analysis, particle swarm optimization, and ImageJ) to design an ordinary differential equation model of caspase-mediated prostate cancer apoptosis signaling. We apply the model to identify key modalities that increase the propensity toward apoptosis across three separate pro-apoptotic drugs (Tocopheryloxybutyrate, Narciclasine, and Celecoxib). Overall, we demonstrate that apoptosis dynamics can be accurately captured in response to each of the three drugs and identify which features of the model represent viable targets for overcoming intrinsic drug resistance.

AbstractsO_ST_ABSBackgroundC_ST_ABSQuantitative systems pharmacology (QSP) models provide mechanistic insight into drug response but are limited by labor-intensive, expert-driven workflows. We developed an AI-assisted QSP (AI-QSP) framework that integrates large language models (LLMs) with SBML-based modeling to enable automated reconstruction, extension, and calibration of mechanistic models. MethodsThe framework was applied to a published CAR-T QSP model. The model was reconstructed in SBML and extended via LLM-guided prompts to incorporate key resistance mechanisms: T-cell exhaustion, PD-1/PD-L1 checkpoint regulation, and tumor antigen escape. Model development followed an iterative expert-in-the-loop workflow. The resulting model (21 reactions, 9 species) was calibrated to synthetic benchmark data using 19-parameter optimization. Model credibility was assessed using ASME V&V 40 and ICH M15 principles, including global sensitivity and profile-likelihood analyses. ResultsThe calibrated model reproduced benchmark dynamics with high accuracy (mean log-RMSE = 0.132). Sensitivity analysis identified CAR-T killing and bystander cytotoxicity as dominant drivers of tumor response. Profile-likelihood analysis showed 71% of parameters were practically identifiable, with remaining parameters prioritised for future data-driven refinement. ConclusionsAI-assisted QSP modeling enables reproducible, scalable model reconstruction and evolution while maintaining mechanistic transparency and regulatory alignment. This framework provides a foundation for accelerating model-informed drug development in cell and gene therapies.

Circadian rhythms are self-sustained biological oscillations that coordinate diverse physiological processes in plants, including growth, metabolism, and environmental responses. These rhythms arise from an interconnected transcriptional translational feedback network that integrates multiple entrainment cues such as light and temperature. The plant circadian clock is organized around key regulatory loops involving CCA1, LHY, PRRs, TOC1, ELF4, LUX, and other transcriptional regulators, whose coordinated interactions ensure precise and robust oscillations. In this study, we developed an ordinary differential equation based mathematical model, building upon a previous framework to incorporate additional regulatory modules and transcriptional controls that better reflect experimentally observed behaviour. To elucidate the regulatory organization of this model, we performed a multi-layered computational analysis combining four complementary approaches: (i) period sensitivity analysis to quantify how parameter perturbations influence the systems timing, (ii) phase portrait analysis to visualize dynamic interactions among key components, (iii) knockout analysis to identify parameters essential for sustained rhythmicity, and (iv) network impact analysis using composite weighted network indices to evaluate hierarchical control across the network. Together, these analyses reveal that transcriptional repression, protein degradation, and light-regulated synthesis form the dominant control mechanisms within the circadian system. The results highlight a hierarchical and robust network structure centred on the CCA1/LHY and PRRs feedback loop, with redundant modules ensuring stability under perturbations. Thus, this model provides an improved, biologically consistent framework for dissecting the dynamic architecture of the plant circadian clock and guiding future experimental validation.

Acquired Dasatinib Resistance (ADR) hinders efficacious treatment of Pancreatic Cancer (PC), often mediated by dynamic signalling reprogramming due to prolonged drug intake. With a novel signalling cross-talk network modelling, this study analyses transcriptomics data of dasatinib-resistant and dasatinib-sensitive pancreatic cancer cell lines and prioritizes key signalling molecules via systemic coordination of their magnitude of dysregulation and the degree of signalling cross-talk among enriched pathways. Results found the p53 and FC-{epsilon} RI signalling pathways demonstrating significant perturbation enrichment, complementarily orchestrating a total of 87% of the global perturbation map in dasatinib resistance. Further statistical characterization of the cross-talk network identified 10 key resistant biomarkers, including THBS1, CDKN1A, and BCL2L1 within p53 signalling, and RAC2 and MAPK13 within FC epsilon RI signalling. Validation with TCGA transcriptomics, CPTAC relative proteomics, and StringDB protein-protein interaction data for their potential prognostics revealed BCL2L1 as pivotal for global perturbation dissemination and, thereby, a novel therapeutic target.

Alzheimers disease (AD) is a multifactorial neurodegenerative disorder characterized by coordinated dysregulation of multiple genes, requiring system-level approaches to elucidate underlying molecular mechanisms. While existing computational studies largely focus on differential expression analysis or machine-learning-based feature selection, they often overlook inter-gene relationships and network topology, limiting biological interpretability. In this study, we present a network-based framework for prioritizing candidate genes in Alzheimers disease by integrating gene co-expression network analysis with multiple centrality measures. Transcriptomic data comprising approximately 39,000 genes across 324 Alzheimers and control samples were preprocessed using log-transformation, variance filtering and Z-score normalization, followed by LASSO-based feature selection to retain phenotype-associated genes. A weighted gene co-expression network was constructed using Pearson correlation to capture coordinated transcriptional activity. Network topology was analyzed using degree, betweenness and eigenvector centrality to identify genes that are highly connected, act as information brokers or occupy influential positions within the network. A consensus ranking was derived by merging these centrality measures, enabling robust prioritization of candidate genes. The results highlight a subset of highly central genes, including several small nucleolar RNAs and regulatory transcripts implicated in RNA processing, synaptic function and neurodegenerative pathways. By jointly leveraging co-expression structure and complementary centrality metrics, the proposed framework provides an interpretable and reproducible strategy for identifying biologically meaningful gene candidates, offering potential utility for biomarker discovery and therapeutic target prioritization in Alzheimers disease.

Recent studies have established that the circadian clock influences onset, progression and therapeutic outcomes in a number of chronic conditions including cardio-metabolic diseases and cancer. For the latter, they also suggest that chronotherapy offers the potential to refine current treatments and improve the development of future anti-cancer medicines. Therefore, there is a need for tools to measure the functional state of the tumoural circadian clock in patients. We have previously developed a model-led machine-learning algorithm called TimeTeller which has the potential to provide such a tool. Here we demonstrate its potential for clinical relevance by a study of 1286 breast cancer patients in which we characterise the nature of the circadian clock disruption in their tumours and demonstrate a strong nonlinear association between 10-year survival and TimeTellers tumoural clock disfunction score {Theta}. We find that good tumour clock function is antagonistic to survival.

The human cell cycle is a highly regulated process that integrates multiple signaling pathways and checkpoints to ensure faithful genome duplication and cell division. Disruptions in these regulatory networks contribute to a wide range of diseases. Here, we present a novel, updateable computational model of the full human cell cycle that shows sustained oscillations over time and reproduces experimental perturbations. We used a hybrid framework combining mass action and Michaelis-Menten kinetics, incorporating the synthesis, degradation, and regulation of key cell cycle proteins and protein complexes. It consists of 63 distinct biochemical species, interacting through 41 major reactions, and functioning through 63 ODEs. The model is built upon a modular framework, structured around the core regulatory networks of the G1, S, G2, and M phases. Due to its complexity, we determined parameter sets that met strict criteria, namely event timing, comparable concentrations, and continuous cycling. We validated the models behavior by reproducing canonical checkpoint responses, including mitogen dependence and the DNA damage response, both of which produced reversible and robust cell cycle arrests. Importantly, the model was trained and calibrated using in vitro data from human U251-MG glioma cells expressing the FastFUCCI cell cycle reporter. We quantitatively aligned the simulated and experimentally determined phase durations and cell doubling times. Next, we experimentally tested and refined model parameters by using abemaciclib-mediated inhibition of CDK4 and volasertib-mediated inhibition of PLK1. In vitro and in silico data show dose-dependent G1 arrest by abemaciclib and dose-dependent mitotic arrest by volasertib. Finally, we demonstrated that the model predicts changes in cell proliferation over a wide range of drug concentrations and combinations. Overall, our work establishes a robust, data-driven computational model for systems-level analysis of the human cell cycle and its disruption by therapeutic perturbations. AUTHOR SUMMARYKnowledge of the protein-protein interaction networks governing the cell cycle is ever-expanding, yet this information is often fragmented across studies focusing on disconnected subsets of the cycle. For decades, researchers have investigated the underlying mechanisms of cell division, but an integrated, quantitative understanding of the entire process remains elusive. This gap is a major hurdle for predicting how targeted therapies affect cell proliferation, especially when used in combination. Our goal is to develop an in silico simulation of the complete human cell cycle by integrating the key mechanistic relationships across all four phases into a single computational model with enough resolution to approximate outcomes upon perturbation. In achieving this, we have developed a novel, comprehensive computational model that provides an integrated quantitative understanding of how cancer drugs alter the human cell cycle. We have rigorously trained, calibrated, and validated the model by suitably estimating its parameter values to produce accurate cell cycle phase timing, using high-resolution, live-cell imaging data and other cardinal features of the cell cycle in a U251-MG glioblastoma line. This work provides an accessible tool for exploring how normal cell cycle control is disrupted in disease, generating new hypotheses, and identifying potential points of therapeutic intervention.

Coral, as a bioreactor, has to continuously interact with surrounding environment to maintain a healthy state. A multi-physics reaction engineering model has been developed to capture this interaction. The coral interior is modeled as interconnected reaction units respectively for photosynthesis, respiration, and calcification, whose reaction kinetics are influenced by environmental fluctuations. Coupling between coral and environment is realized by bi-directional mass transfer at the coral-seawater interface, with consideration of the unique flow fields induced by ciliary beating. By resorting to this comprehensive model, we discover that ciliary beating demonstrates distinctively different diurnal and nocturnal functions. During daytime, beating can help reduce photosynthetic oxygen accumulation to prevent hyperoxia-induced mortality, while enhancing carbon dioxide uptake efficiency to promote nutrient production. At night, however, beating promotes oxygen acquisition for adequate respiration, while expelling carbon dioxide to inhibit symbiotic destruction under acidic stress. The model further enables mechanistic analysis of the detrimental impact of climate change on coral health, where the influences from two key factors (i.e., temperature and CO2 level) can be decoupled. Its interesting to find out that the elevated temperature plays a dominant role during daytime, while at night the coral is dominantly influenced by rising CO2 level.

BackgroundFor decades, computational biology has failed to create unified models where metabolic state and regulatory control are bidirectionally coupled: metabolic models optimize flux but cannot represent dynamic regulation, while regulatory models treat ATP as a fixed parameter rather than a dynamic variable affected by pathway activity. This fundamental limitation prevents computational recapitulation of emergent threshold behaviors--spontaneous homeostasis, adaptive reorganization, pathway switching--observed in living organisms. The challenge requires formalisms where (1) metabolic state governs regulatory decisions AND (2) regulatory choices consume metabolic resources, producing emergent dynamics from feedback rather than programming. MethodsWe introduce Signal Hierarchical Petri Nets, extending Hybrid Petri Nets with bidirectional metabolic-regulatory coupling through energy-dependent layer organization. Unlike classical approaches, ATP is simultaneously a regulatory signal (governing pathway availability through quantitative thresholds) and a material substrate (consumed by pathway activity). When ATP depletes below 1000 {micro}M, high-cost pathways automatically become unavailable; pathway activity consuming ATP creates feedback affecting subsequent pathway accessibility. This bidirectional coupling enables emergent threshold behaviors impossible in classical formalisms. We demonstrate the paradigm through macrocyclic peptide transport across 53 metabolic conditions, where drug accumulation depends on ATP-governed pathway reorganization. ResultsThe formalism produces three emergent behaviors never achieved in unified metabolic-regulatory models. (1) Spontaneous homeostasis without programming: Despite 113-fold permeability variation from N-methylation, ATP-replete cells maintain constant drug accumulation (CV=0.066%)--homeostatic compensation emerges from ATP-consumption feedback, not explicit control logic. (2) Threshold-triggered reorganization: ATP depletion to 300 {micro}M triggers 8533-fold active-to-passive transport shifts with paradoxical 141% accumulation increase from efflux collapse. (3) Tissue-specific dynamics from identical parameters: Tumor hypoxia (ATP=1200 {micro}M) versus normal tissue (ATP=5000 {micro}M) produces 6.62-fold selectivity differences from differential pathway accessibility--same model, different emergent outcomes. Computational predictions achieve r=0.911 correlation with experimental cyclosporin permeability (n=32). ConclusionsSignal Hierarchical Petri Nets represent the first computational formalism achieving emergent threshold dynamics through bidirectional metabolic-regulatory coupling. The paradigm enables in silico recapitulation of adaptive cellular behaviors previously impossible to model, with applications extending beyond drug transport to any biological system where metabolic state governs regulatory reorganization: cancer metabolism, ischemia, synthetic biology, and aging research. Author SummaryLiving cells exhibit remarkable adaptive behaviors: they maintain stable internal conditions despite environmental changes (homeostasis), reorganize their biochemical machinery when energy runs low, and switch between operating modes at precise threshold values. For decades, computational biologists have struggled to build models that recapitulate these emergent behaviors--our computer simulations could only exhibit the dynamics we explicitly programmed into them. We solved this fundamental challenge by creating a new computational formalism where metabolic state and regulatory control form a bidirectional feedback loop: energy availability governs which biochemical pathways can operate, while pathway activity consumes energy. This simple coupling produces complex emergent behaviors never seen in previous computational models. Our simulations spontaneously exhibit homeostasis--maintaining constant drug levels despite 113-fold variation in membrane permeability--without any programmed control logic. The same model produces different emergent behaviors in different tissue contexts: tumor cells versus normal cells show 6-fold differences in drug accumulation from identical parameters, purely from different starting energy levels. We demonstrate this paradigm using drug transport as a case study, but the implications extend far beyond: cancer metabolism, brain injury during stroke, synthetic biology circuit design, and aging research all involve systems where metabolic-regulatory feedback governs cellular adaptation. This formalism provides computational biology with a long-missing capability--the ability to model emergent threshold behaviors computationally.

Drug discovery is impeded by the difficulty of translating targets from preclinical models to patients. In this work, we present TwinCell, a Large Causal Cell Model for target identification that, trained on in vitro cancer cell line perturbation data, generalises to patient-derived cell types while providing biologically meaningful interpretations of its predictions. Rather than predicting perturbation outcomes, TwinCell identifies the upstream regulators most likely to drive the transition between two cell states, such as diseased and healthy, by decomposing target probability over signalling paths through a multiomics interactome conditioned on single-cell foundation model embeddings. To validate this approach, we introduce TwinBench, a benchmarking framework that evaluates virtual cell models using recommendation system metrics while correcting for mode collapse through empirical p-value estimation. On both in vitro zero-shot scenarios and in clinico validation across five therapeutic areas, TwinCell outperforms not only state-of-the-art virtual cell models but also linear baselines and network-based methods, classically used to perform target identification. When applied to patient data, TwinCell recovers clinically approved drug targets and reconstructs known disease mechanisms, such as the type I interferon signalling cascade in Systemic Lupus Erythematosus, without any disease-specific supervision. These results demonstrate that constraining learned perturbation patterns to a biological interactome enables cross-tissue, cross-disease target identification with mechanistic interpretability, bridging the gap between high-throughput in vitro experiments and clinical insights.

Metabolic reprogramming is a hallmark of glioblastoma, yet how distinct malignant and tumor microenvironment cell populations contribute to this metabolic heterogeneity remains poorly defined. Since direct single-cell metabolomics remains technically limited, transcriptomics-based computational inference offers a powerful alternative. Here we apply and systematically compare three complementary computational methods: (1) metabolic pathway activity scoring, (2) gene regulatory network inference focused on metabolic enzyme gene regulation, and (3) single-cell metabolic flux prediction. These methods were applied to snRNA-seq data from a set of GBM patient samples using the Human1 genome-scale metabolic model as a unified reaction and pathway annotation prior knowledge reference. Across all three methods, tumor-associated macrophages emerge as the metabolically dominant tumor microenvironment population. Tumor-associated macrophages in mesenchymal-like tumors show coordinated transcriptional control of lipid metabolism by five recurrently active transcription factors. They also exhibit consistent nucleotide biosynthesis flux and glutamate-to-glutamine conversion potentially supporting malignant cells. These findings demonstrate that multi-layered metabolic inference can resolve cell-type/state-specific dependencies in glioblastoma and highlight tumor-associated macrophage metabolism as a promising therapeutic target

The pancreas is a complex inner organ that develops through the concurrent formation of multiple interacting biological networks, but the geometric principles governing their spatial organization remain poorly understood. Here, we study pancreatic development by analysing loop structures in ductal, neuronal, and vascular networks using tools from topological data analysis. In particular, we employ the application of chromatic persistence, which is designed to detect topological interactions between distinct subsets of a structure; this is used to quantify how loops from different networks become spatially entangled. Our analysis reveals distinct developmental timelines for loop formation, progressive entanglement between ducts and vasculature, and spatial concentration of threading toward the organ interior, providing a geometric characterisation of coordinated network development in the pancreas.

BackgroundThe epidermal growth factor receptor (EGFR) orchestrates highly context-dependent intracellular signaling networks whose architecture varies across cell types and is frequently rewired by targeted therapeutics. Systems-level reconstruction of these networks from phosphoproteomic data remains challenging because phosphorylation measurements identify signaling nodes but do not reveal the interaction paths that propagate signals between proteins. ResultsWe developed a computational framework integrating time-resolved phosphoproteomics with graph traversal algorithms to reconstruct EGFR-initiated signaling pathways across three contexts/conditions. A sign-assignment preprocessing procedure converts quantitative phosphorylation measurements into binary activation states across time points, defining a condition-specific active node set that filters the protein-protein interaction network. Breadth-First Search combined with interaction-weighted Beam Search is then applied to the STRING interaction database (v11.5) to enumerate candidate signaling paths. Applying this framework to phosphoproteomic datasets from EGF-stimulated HeLa cells, EGF-stimulated MDA-MB-468 triple-negative breast cancer (TNBC) cells, and EGF-stimulated MDA-MB-468 cells pretreated with the SHP2 inhibitor SHP099 yielded 260 paths in HeLa cells (117 unique topologies), 293 paths in MDA-MB-468 cells (155 unique), and 292 paths under SHP2 inhibition (85 unique). HeLa cells displayed a SRC-centered architecture dominated by ERBB2 and SHC1 first-hop effectors, converging on focal adhesion, HSP90 chaperone, CRKL adaptor, and integrin signaling arms. In contrast, MDA-MB-468 cells showed a PIK3CA/PTPN11 dual-axis architecture integrating direct PI3K engagement with SHP2-mediated GRB2-IRS1-ABL1 signaling. SHP2 inhibition abolished PTPN11-mediated pathways and induced PIK3CA dominance (69.2% first-hop), accompanied by compensatory ERBB3 engagement and a computationally predicted SYK/VAV1/LCP2 node set whose biological role warrants experimental validation. ConclusionsTime-resolved phosphoproteomics-guided BFS Beam Search over STRING interaction networks captures cell-type-specific EGFR signaling architectures and drug-induced pathway rewiring. This framework provides a systematic approach for transforming phosphoproteomic measurements into mechanistically interpretable signaling hypotheses specific to the cell-type-specific contexts, directly applicable to drug resistance modeling.

Difficult-to-treat (D2T) rheumatic disease affects approximately 12% of rheumatoid arthritis patients and resists sequential biologic therapy, yet no mechanistic model explains this refractoriness as a system-level phenomenon. Here we present the 3-Axis Integrative Framework (3-AIF), a six-variable ordinary differential equation system integrating mucosal tolerance, energy-gated neuroimmune danger sensing, and integrated stress response pathways coupled through Hill-function metabolic gating. Stability analysis reveals bistable dynamics with two co-existing attractors separated by a saddle point. Bifurcation analysis demonstrates fold catastrophe with hysteresis: recovery requires greater therapeutic effort than disease prevention. Sensitivity analysis identifies three dominant parameters mapping to neuroimmune activation, energy drain, and recovery capacity. Cross-disease transcriptomic consistency analysis across six public datasets (n=310, five rheumatic diseases, four tissue types) reveals compartment-specific axis dysregulation -- circulating cells show integrated stress response activation while target tissues show pathway exhaustion -- and disease-specific axis dominance patterns consistent with model predictions.

Spatial computational models such as agent-based models (ABMs) offer powerful in silico tools to study tumor dynamics, yet imaging data are still rarely used to inform these models directly. We present an ABM optimization framework that leverages convolutional encoders to compare spatial patterns between experimental imaging data and ABM-generated outputs within a shared latent space. This quantitative comparison was used to estimate ABM parameters across three datasets, ranging from synthetic data to 3D tumoroid-T cell co-culture microscopy and histopathology images from The Cancer Genome Atlas skin cutaneous melanoma samples. Estimated parameters were evaluated using data-derived features and experimental knowledge, including experimental conditions and gene expressions. Simulations using optimized parameters reproduced key spatial features of the training images, such as tumor boundary complexity and tumor-tumor neighborhood structure. Together, these results demonstrate a flexible framework for ABM parameter optimization using spatial data across modalities, enabling systematic investigation of how spatial architecture influences tumor progression and immune interactions.

Agent-based modeling (ABM) is a computational method for predicting the emergent outcomes of interacting, autonomous individuals in a complex system. Here, ABM is used to simulate interactions between fibroblast and myofibroblast cells during idiopathic pulmonary fibrosis (IPF) in alveolar tissue microenvironments. These microenvironments are derived from histology of a healthy human lung sample and moderate- and severe-IPF lung samples. Fibroblast differentiation, cell migration, and collagen secretion in response to the spatial distribution of the cytokine transforming growth factor-beta are captured in the ABM using NetLogo software. Results are presented from one simulated year without treatment and with mechanisms representing treatment by pirfenidone and pentoxifylline, alone and in combination. A total of 180 in silico experiments are run, analyzed, and compared in a high-throughput workflow. The effects of the initial number of fibroblasts and treatment scenarios on various metrics related to collagen accumulation and collagen invasion into alveolar regions are determined. The ABM and the analysis files are shared to facilitate model reuse. By integrating computational modeling of IPF and therapeutics, this research aims to improve understanding of fibrosis progression and assess the efficacy of novel and existing treatments targeting different mechanisms to inform decision-making for IPF treatment.

Understanding how gut microbes transform drugs, and how this influences microbiome composition and function, remains a key question to better understand the efficacy and side effects of pharmaceuticals. To accelerate the discovery of microbiome-derived drug metabolites, we developed a functional metabolomics workflow that combines the use of synthetic microbial communities (SynComs) with a time-series resolved molecular networking approach and advanced computational metabolite annotation. We demonstrate how this framework can be used to illuminate chemical transformation dynamics in a gut SynCom (Com20) with 50 clinical drugs. Our results highlight a multitude of drug metabolites, including multi-step metabolic cascades, some of which correlated to shifts in microbial taxa, suggesting functional links between microbiome composition and biochemical transformations. Our computational data analysis workflow is publicly available through the GNPS2 ecosystem at chemprop.gnps2.org, which can be used to prioritize biotransformations and other (bio)chemical reactions in various biological and abiotic systems.

Dysregulated kinase activity drives oncogenic signaling, disrupts cellular homeostasis, and promotes tumour progression. The BRAFV600E mutation constitutively activates the MAPK pathway and is a key therapeutic target in melanoma and other cancers, but the functional relevance of most downstream phosphorylation events and mechanisms of drug resistance remain unclear. To address this, a global multi-omic model of BRAF inhibition response was established in BRAFV600E-mutant cells by integrating time-resolved and biophysical phosphoproteomics, transcriptomics, and thermal proteome profiling. Ultradeep phosphoproteomics revealed extensive phosphorylation changes upon BRAF inhibitor treatment, while biophysical phosphoproteomics identified phosphorylation events linked to altered protein solubility and subcellular localization, suggesting changes in nucleic acid interactions and nuclear reorganisation. Network-based integration of these datasets prioritized functionally relevant phosphorylation sites and kinases. Experimental validation identified CDK9, CLK3, and TNIK as critical regulators of BRAFV600E signaling and candidate targets for combinatorial inhibition capable of re-sensitising resistant cells. The transcription factor ETV3 emerged as a previously unrecognised effector of BRAF signaling. Biophysical proteomics data confirmed that ETV3 phosphorylation modulates DNA-binding, while functional assays combining knockdown, metabolomics, and drug screening demonstrated its role in coordinating transcriptional and metabolic adaptations to BRAF inhibition. This study provides a systems-level framework linking phosphorylation dynamics to protein function and phenotype, identifies ETV3 as a new node in oncogenic BRAF signaling, and illustrates how integrated, site-resolved models can reveal mechanisms of kinase-driven oncogenesis. O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=43 SRC="FIGDIR/small/704793v1_ufig1.gif" ALT="Figure 1"> View larger version (13K): org.highwire.dtl.DTLVardef@1df1401org.highwire.dtl.DTLVardef@9a77a5org.highwire.dtl.DTLVardef@124f819org.highwire.dtl.DTLVardef@1c6b57_HPS_FORMAT_FIGEXP M_FIG C_FIG HighlightsO_LITime- and cell type-resolved phosphoproteomics maps BRAF inhibition dynamics C_LIO_LIBiophysical phosphoproteomics, combining quantitative phosphoproteomics with solubility profiling or nuclear fractionation, reveals phosphorylation-driven changes of protein solubility and localization C_LIO_LIIntegration of abundance and biophysical phosphoproteomics data identifies functionally relevant phosphorylation events of BRAFV600E signaling C_LIO_LINetwork integration of multimodal phosphoproteomic, transcriptomic and thermal proteome profiling data links signaling to protein function and cellular phenotypes C_LIO_LIBiophysical evidence improves models and identifies non-canonical kinases driving BRAF signaling as well as novel downstream regulators such as ETV3 C_LIO_LIFollow-up experiments reveal a ETV3-GLUT3-mediated metabolic adaptation in BRAFV600E cells C_LI

Stachys affinis (Chinese artichoke) tubers contain 50-80% stachyose by dry weight, the most concentrated dietary source of raffinose-family oligosaccharides (RFOs) known. Because humans lack sufficient -galactosidase activity, stachyose transits intact to the colon where microbial fermentation yields short-chain fatty acids (SCFAs). However, the quantitative impact of stachyose-derived SCFAs on host hepatic energy metabolism has not been systematically explored using genome-scale metabolic models. Three stachyose dose scenarios (Low/Mid/High: [~]25, 50, 100 g fresh tubers) were translated to SCFA availability vectors. Hepatic metabolic responses were simulated using Recon3D (10,600 reactions) and Human-GEM (13,417 reactions) under strict hepatocyte-like media, maximizing ATP maintenance flux (ATPM). FVA across multiple optimality thresholds (90-100%) and pFBA confirmed solution robustness. One-at-a-time sensitivity analysis characterized ATPM responses to individual parameter perturbations, and a ratio sensitivity sweep across six alternative SCFA profiles assessed dependence on assumed fermentation ratios. A targeted rescue experiment addressed model-specific propionate catabolism gaps. Both models showed dose-dependent ATPM increases (Recon3D: +71 to +286%; Human-GEM: +103 to +413% above baseline), with the 19-33% inter-model gap attributable entirely to Human-GEMs functional propionate catabolism pathway. FVA confirmed near-unique optimal solutions (ATPM ranges [~]1% at 99% optimality, widening to [~]10% at 90%). Parsimonious FBA preserved identical ATPM values while reducing total flux by [~]4-14%, confirming objective robustness. SCFA ratio sensitivity across six alternative profiles showed 27- 28% ATPM variation, indicating qualitative robustness. Butyrate yielded the highest ATP per mole ([~]22) in both models; propionate sensitivity was zero in Recon3D but [~]15.25 mmol ATPM/mmol propionate in Human-GEM. Reopening propionyl-CoA carboxylase (PPCOACm) in Recon3D under strict constraints converged ATPM to within 0.3-0.7% of Human-GEM, cross-validating both reconstructions. This reproducible dual-model pipeline identifies model-specific pathway gaps and provides cross-validated predictions to guide future experimental studies of how dietary SCFAs influence hepatic ATP metabolism.